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Abstract 



We propose a model describing the liquid-vapour phase transition according to 
a phase-field approach. The model takes up the setting proposed in [1], where 
O^l a phase field (p is introduced whose equilibrium values ip — and ip — 1 are 

' associated to the liquid and vapour phases. The phase field obeys Ginzburg- 

Landau equation and enters the constitutive relation of the density, accounting 
for the sudden density jump occurring at the phase transition. In this paper 
we concern ourselves with the extension of this model to take into account the 
fH existence of the critical point in the coexistence line. 

> 

^2 Keywords: Liquid- Vapour Phase Transitions, Critical point. 

Thermodynamics. 
I 1 PACS: 44.35. -l-c, 64.70.fm, 68.35.Rh. 

0^ 1. Introduction 

\o 

In this paper an idea put forward in [T] concerning a phase field description 
of water-vapour phase transition is developed. The main macroscopic manifes- 
04 tation of the liquid vapour phase transition (and, to a smaller extent, of the 

ice-water transition) is a sudden jump at a given pressure in the density of the 
substance. The starting assumption made in the cited paper is the decomposi- 
tion of the density p of the fluid into two contributions po and pi, such that 



■>< i = ^ + ^. (1) 



The contribution po is taken as a function po = PqIPt^) of the pressure p and 
temperature 9, while pi also hinges on an independent phase field ip. Such a de- 
composition entails the attainment of the Andrews curves with their character- 
istic plateaus. Otherwise, these would be obtained from the van der Waals state 
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equation supplemented by the Maxwell construction. In the cited approach, in- 
stead, it is the variation of the order parameter which gives rise to the horizontal 
piece of the isothermal lines (which prevents the relation between pressure and 
specific volume to be one-to-one). The additional phase field requires a further 
evolution equation (besides continuity, Navier-Stokes and heat equations) which 
is of course the time-dependent Ginzburg-Landau equation. We remark that in 
the classical Ginzburg-Landau approach to the liquid-vapour critical point, the 
order parameter is identified with the difference p — pc between the density p and 
the critical density pc (see e.g. 2, 3]). Here, the order parameter is independent 
of the density of the fluid and stands out from the constitutive decomposition 
([I]) of the overall fluid density. 

In this paper we concern ourselves with the modeling of the liquid-vapour phase 
transition along the same lines of [1], but also including the peculiarities due 
to the presence of the critical point. This aspect was not fully considered in 
the previous model, which makes sense specifically for temperatures below the 
critical temperature. As a consequence of the existence of the critical point, 
there always exists a sequence of equilibrium states which connect in a contin- 
uous way the liquid and the vapour phase (see Fig. [T]). Taking into accoimt 
this kind of transformations, it is not possible to assign a unique fixed value 
of ip to the liquid and vapour phases anymore. We remark that, due to the 
existence of the critical point, it is possible to speak of distinct phases only on 
the coexistence curve. So it is necessary to consider a phase field which admits 
a continuous set of equilibrium values under different temperature and pressure 
conditions. Obviously this exhibits the second order (or continuous) character 
of the transition at the critical point. In usual phase field models (see e.g. [Ull]) 
the homogeneous phases of a substance are associated to fixed discrete values 
of the order parameter (for example if = and if — 1). Intermediate values 
are associated to the thin interfacial layers between two phases; in particular, 
they are associated, in stationary conditions, with non vanishing values of \7ip. 
This will not be true anymore if we assume a phase field which have a contin- 
uous set of equilibrium values associated to homogeneous phases, similarly to 
a Ginzburg-Landau order parameter in a strict sense (for second order phase 
transition) . In this case intermediate values of the phase field can be associated 
both to interfacial layers and to homogeneous phases, close to the critical point. 

2. Model equations and thermodynamical consistence 

We start by setting the general thermodynamic framework of the model we 
are going to propose. Our description is based on the introduction of an order 
parameter (p connected with the density jump of the transition. The equation 
for if is formulated as an additional balance equation to which is associated a 
corresponding balance of powers [J [51 15] . As a consequence, a variation of pi is 
associated to a supplementary power flux in addition to the classical mechanical 
power pp/p. 

We construct an evolutive model based on four balance equations for the 
fundamental thermomechanical flelds of the fluid, that is: the phase field ip, the 
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density p, the macroscopic velocity v and the temperature 9. 

The thermomechanical state of the system is defined by the set of variables 



a= (p,0,<^,h), (2) 

where h = V</?. We remark that we are considering a compressible fluid, both in 
its liquid and in its vapour phase. In this situation, the choices of the pressure or 
the density p as state variables are equivalent. From our point of view, it is more 
natural to choose the pressure as state variable, because it appears explicitly in 
the Ginzburg-Landau equation for the order parameter. 

We state the four balance equations of the model in the order: order balance 
equation, mass balance, momentum balance and energy balance 

p/c - V • p, (3) 

p + pV • V = 0, (4) 

pv = Vp + V • Tb +pb, (5) 

pe^V' + [}, (6) 



where: 



the superimposed dot represents the material derivative of a field: (j) = 
d(j)/dt + V • V0; 

k is the internal order production and p the order flux. For the meaning 
of this terminology we refer to [S]. The relation of these flelds with the 
order parameter ip is given in Eq. ([9]); 



T^; is the extra stress (see Eq. (14|), and b the external force density; 



e is the internal energy density, 7" is the internal power density which 
is the sum of the internal mechanical power associated to motion 
equation and the power T'^ associated to order balance; 

f) is the internal thermal power density, which obeys its own balance in 
terms of heat flux q and heat supply r: 

f)^-V-q + r. (7) 



The balance laws have to be completed by appropriate constitutive relations for 

^,P,T£;,p. 

The phase field equation suitable for the description of a phase transition 
takes the form of a Ginzburg-Landau equation: 

pTip = kV • (pV^) - pf^{p, 0, ip) (8) 

where r and k are positive parameters, / is the Ginzburg-Landau potential and 
ftp = df /dip. Hence the fields k and p satisfy the constitutive equations 

k^Ttp + f^, p = KpVtp. (9) 
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Equation ([8]) allows us to write the balance of powers associated to order 
parameter by multiplying both members for ip and extracting a divergence con- 
tribution: 

pTip^ + pfy,(p + p\7(p ■ \7{Kip) = V • {npifS/ip). (10) 

We identify the first member with the internal power of the order parameter; 
using the identity 

V(^ = h + h-Vv. (11) 

we find useful to write the internal power as 

Vl, = p{Tp^ + UP + Kh-h + Kh-Bh), (12) 

where D = ^(Vv + Vv-^) is the rate of deformation tensor. 

The internal mechanical power V^, has the classical expression 

n.-/+TB:D (13) 

and the extra-stress is assumed in the form 

Te = ^i{^,e)B - Kph(»h. (14) 

It includes a viscous contribution and a coupling term with the order parameter 
gradient h which is necessary to assure the thermodynamical consistence (see 
below). The viscosity fi{6,(p) is a function of the temperature and the phase 
which vanishes in the vapour and gas phases. 

The constitutive equation for p is conveniently expressed in terms of the 
specific volume v := 1/p. We will assume a generic constitutive relation v = 
uip, 6*, (fi) The choice of the function i> is strictly related to the function 
f{p,9, (p) by the thermodynamical constraints and will be specified after we will 
have examined the consequences of the Second Law. We find it convenient to 



express the mechanical power (13) in terms of the specific volume i> 



= -pPii'pP + i>eO + v^ifi) + Te :T>. (15) 

Now we examine the constraints on the constitutive equations due to the 
Second Law, in the form of the Clausius-Duhem inequality. In terms of the free 
energy ip = e — Orj, the inequality can be expressed as 

p{^ + TjO) - - + • g < 0, (16) 



where g = V9. Substituting the expressions ( 12 ) and ( 15 ) we obtain 

(V'e + 77 + pi>0)d - Tip^ + {ijj^ - U+ pi'v)^ + (V'p + pi'p)? + (17) 

{iph - Kh)h - {-Te + Kh ® h) : D + • g < 0. 
P P0 



^This is different from the constitutive relation for density, pressure and temperature which 
is assumed in which is rate dependent. 
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Firstly we observe that constitutive equation ( 14 ) makes the term proportional 
to D a non- negative definite term, provided that fj- > 0. From inequality (17), 
it is plain the reason to include the contribution proportional to h (g) h in the 
extra stress tensor. 

In order to satisfy the Clausius-Duhem inequality in the most simple (non 
trivial) way, we add the following constitutive relations 

1] = -ipe-pug, (18) 

V'lp = f^-pi'v, (19) 

pi>p = -il)p, (20) 

V-h = Kh, (21) 

q = -fcog, fco > 0. (22) 

Then, the following expressions are easily obtained]^ 

V;(p,0,¥',h) = ^h2 + /(p,0,^)-p/p(p,0,^), (23) 

v{p,e,^) = fp{p,e,ip), (24) 

viP,0,ip) = ~feip,0,ip). (25) 

Note that 

$ ^h^ + / (26) 

takes the meaning of a specific Gibbs' free energy and ■(/; = $ — p^p. 

Substituting the expressions of D and in the internal power expression 
we obtain 

■p* - p(V + + + /i||D|p. (27) 
Then, the first law of thermodynamics 

pe = -V • q + pr + -p* (28) 

gives us the heat equation 

-pOife)- ^V-ikoV9)+pr + pT^^+fx\\T>\\^. (29) 



We resume the set of equations describing the model 

prip ^ K-W{pW(p) - pf^{p,9,(p) 
/9 + pV • V = 

< = Vp + V • {fi{(p, 6')D - Kph <^h)+ ph (30) 
-p0{fey - V • (kaVe) + pr + prip^ + ^||D||2 



■^Actually, one would obtain, for a given /, that VCPi Vi h) = Kh^/2 + /(p, 6, !p)+g{p, 8) — 
pfp, but, as g does not depend upon ip, we could rename / + g as /, without modifying 
Ginzburg-Landau equation. 
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The last equation is the constitutive equation for the density. It could be re- 
garded as constitutive equation for the pressure if we would consider the density 
as a state variable. 

The Clausius-Duhem inequality turns out to be very useful in establishing the 
coupling between the different equations. It has nothing to say about the choice 
of the function / instead. 

Equilibrium properties. 

From the Ginzburg-Landau equation ([s]), the equilibrium value cp = (p{p,6) of 
the order parameter is a solution of the equation 

/^(p,0,^) = O, (31) 

which also should be a local minimum of /, for stability reasons. As the system 
has at most two coexisting phases, f{p,9,ip) has to admit, as a function of ip, 
at most two local minima (^i, (p2- On the coexistence line 

p = pa{e), e<ec, (32) 

the two local minima have to be absolute minima, that is: 

f{ipi,Poi9),9) = f{^2,Po{O),0). (33) 

In general, there is a strip-like region bordering the coexistence line in which 
both minima are present: one of them is an absolute minimum and represents 
the true equilibrium state in that region, while the other one is a relative mini- 
mum and represents a metastable state, which can be a superheated fluid or an 
undercooled vapour. Elsewhere, the order parameter has a unique equilibrium 
value, in particular for 9 > 9c (or p > pc). 



The equilibrium Gibbs free energy ( 26 1 for the homogeneous substance (that is 
for /i = 0) in the z-phase is 

^,{p,9):= f{p,9,^,(p,9)) (34) 

and the coexistence line is characterized by i>i = ^2- 

The usual thermodynamics of phase equilibrium is then recovered. Firstly, we 



remark that, thanks to (31 1, we have 



— g;^ — = fp(P,9,(p)\^=f,,, — 1^ — = fe{p,9,(p)\^=g,,. (35) 



So, the equilibrium equation of state of the i-phase which follows from ( 24 ) is 
represented by the usual relation 

Vi{p,9)^ &p ■ ^ ' 



Similarly, from (25), the equilibrium entropy is given by 



dh{p,9) 

m{P,d) = TTT, • (37) 
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The specific volume jump due to phase change 

Aiy{e) = i>2{po{O),0) - i>i{PoiO),0) (38) 

is one of the most important parameters of the transition, together with the 
latent heat 

L{0) = 0A7^ = e[Upo{O),0)~m{po{9),0)]. (39) 
These quantities are related by the well known Clausius-Clapeyron equation 

L^0p[,{0)Aiy, (40) 

which can be easily obtained deriving both members of equilibrium condition 

^,{poi0),0)^^2ipo{e),e) (41) 



with respect to 9 and using (|36[) , ( 37 1 



3. The Ginzburg-Landau potential 

We now turn to the choice of the Ginzburg-Landau potential /. The phase 
field approach to first order phase transitions involves an order parameter 
which assumes a fixed value (say (p — 0) in one of the two phases, and a dif- 
ferent fixed value (say = 1) in the other phase. The order parameter does 
not suffer jump discontinuities at phase boundaries, but varies smoothly across 
a interfacial layer from one equilibrium value to the other O [7l |8] . This is ob- 
tained by subjecting the order parameter to a Ginzburg-Landau equation, with 
a potential f{ip) having minima only at (/s = 0, 1 for all values of the temper- 
ature and pressure (or any other relevant physical field). In the paper [T] it 
was chosen a potential such that it had always two relative minima at ip = 
and if = 1. That minima were associated respectively to the liquid and vapour 
phases. This setting is not suited to model the phases in the neighborhood of 
the critical temperature. In fact, at the critical temperature 0c, the liquid and 
vapour phases merge in a unique phase which is sometimes named gas phase 
in a strict sense (or supercritical fluid in the region > 0c, p > Pc)- More 
precisely, at the critical point, a second order phase transition occurs. We can- 
not associate a fixed value of the order parameter to the states lying on one 
side of the coexistence line, if this one have an end point. In fact, consider 
the transformation represented in Fig. [T] No phase transition occurs, so the 
order parameter has no jump discontinuity along that way. Therefore the order 
parameter should have a continuous set of equilibrium values. In particular, 
this rules out an approach based on three phases, liquid, vapour and gas, as- 
sociated to fixed values of the order parameter (for example, 1, -1, 0): in that 



We distinguish gas from vapour by saying that vapour can undergo a liquefaction transi- 
tion at sufficiently high pressure, while gas does not. 
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physically inexistent first order phase transition between gas and liquid 
and between gas and vapour would take place somewhere in the phase diagram. 
The most direct generalization of the approach followed in [1] is to assume a 




Figure 1: The phase diagram near the liquid- vapour critical point. The transformation from 
state 1 to state 2 through the path indicated by the arrow does not involve discontinuities of 
any kind 

Ginzburg-Landau potential with minimal points varying continuously with tem- 
perature at 9 < 9c. The two minima which are present in the coexistence region 
of liquid and vapour merge in a unique minimum ip = for 9 > 9^- The pressure 
will not move the position of the minima, but can make the order parameter 
jump from the liquid minimum to the vapour minimum and vice versa. So, in 
isothermal condition and inside the one phase region, the order parameter will 
be effectively constant and the Navier-Stokes equation is in fact decoupled from 
Ginzburg-Landau equation. 

We start by considering the double well potential 

4 2 

F(a:;u) = ^-w2^, u>0. (42) 

having two minima at a; = ±u and a maximum at a; = 0. We next define the 
piecewise polynomial and odd function 

{-^ — u^x if \x\ < u 

(43) 
~^sgn{x)u^ if \x\ > u 

Any linear combination F{x;u) + hG{x;u), with a real coefficient h, is a 
function of x having a minimum or a flcxus in a; = ±u and no other minimum. 
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The Ginzburg-Landau free energy density is taken in the form 



/(p, 0, cp) = Up, 9) + F{p- u{e)) + Hp, e)G{^- u{e)), (44) 

where /i is a function vanishing on the coexistence hne, namely h(po{0), 9) = 0. 
Then, the specific volume is given by 

v{p,9,^) = fp{p,9,^) = {h)p[p,9) + H[p,9)G[^-u{9)). (45) 

We identify the components of the decomposition ([T]) in the following way 

- = (/o)p(p, 0), - = hp{p, 9)G{ip- u(9))- (46) 
Po Pi 

Then, the difference in volume on the coexistence line, corresponding to the 
jumping in Lp from the minimum (pi ~ u io the minimum (p2 = is 

^y^^hM(^)^e)u\ (47) 

If we assume h{p,9) = p — Po{9) for 9 < 9c, then the parameter u is directly 
related to the measurable volume difference 

u{9) = i^-^AH9)j . (48) 

The liquid phase corresponds to the phase diagram region 9 < 9c, p > Po{9), 
that is, in the case we are considering, to /i > 0. For h > the absolute minimum 
is attained a.t ip = u, so, for this choice of h, the liquid phase is associated to 
the right minimum, while the vapour phase is associated to the left minimum. 
For 9 > 9c we have Ai/ = and so u{9) = 0. We note that G{x; 0) = 0, so, for 
9 > 9c the Ginzburg-Landau potential is 

f{p,9,p)^foip,9) + ^ {9>9c), (49) 

which has the only minimum p = Q. We associate that value to the gas phase. 

On isothermal and homogeneous conditions (so that V(/3 = 0) and with ex- 
ternally prescribed homogeneous pressure field p{t), the evolution of the volume 
is ruled by the equations 

T<^ = -{F^{<p;u) + h{p{t),9)G^{p,;u)), (50) 
Ht) = {.fo)pipit),9) + hp{p{t),9)G{p;u). (51) 

For T ^ 1, one can assume that at every instant the order parameter is equal 
to the equilibrium value (p{t) = ±u, so that the equilibrium equation of state is 

Hp, 0) = (/o)p(p, 9) + hpip, 9)G{p; u). (52) 
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h A 




Figure 2: Structure of the minima of the potential f{p,8,ip) at varying h > 0,u > 0. In the 
region h < the structure is the same with the interchange of ip with —ip. 



The equilibrium value (p considered as a function oi p,6, is not a single- valued 
function, since there is an hysteresis behaviour in the interval \h\ < u. Con- 
sidering only the stable states and discarding the metastable ones, amounts to 
take (p as 

r u ifp>po{0) 
^ip,e)^l (53) 
[ ifp<po{9) 

Assuming this function ip{p, 0) the equilibrium isotherm at 9 < 9c in the 
u — p plane is 

i^ip,9)^l (54) 

which has a flat piece at pressure po{9) with length lS.v{9) = Av? /"i. For pres- 
sures p ^ po(^) the curve is determined essentially by /o(p, 9). For 6 > 9c, if = 
and ly = (/o)p. 

A difficulty in this model is given by a possible diverging contribution in 
the equilibrium entropy if Ai'{6) is singular at 6* = 6'c, as, for example, if Az^ ~ 
\9 — 9c\^ , with /3 < 1. In fact, consider the expression of the entropy 

r]{p, 6*, ip) = -fe{p, 9, ip) = r]o{p, 9) + uugip^ - hgG{ip; u) + 2huug(p, (55) 
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where r]o(j),9) :— So, a diverging expression for ue would give rise to 

a diverging entropy at critical point. We can directly relate the equilibrium 
entropy fj with Au{9). At the equilibrium values (p — ±u, we have the entropies 

2 

fj± = rjo + U0U^ ± -hgu^ ± 2hu^ug, (56) 
or, remembering that Ai/((?) = 4u'^/3, 

^±=^o+(^±0(A.)«±^A.. (57) 
In particular, at the critical point u = 0, Av = and the entropy is given by 

fi±ip, Oc) = Vo{p, Oc) ± hi9,,p)iAiy)e. (58) 

Now, in general {Ai')q is a singular quantity at the critical point; for exam- 
ple, according to the Landau theory of the second order phase transitions, 
Ai> cx |0 — So, in this case we would obtain a divergent entropy at 

the critical temperature. In paper ^ both equilibrium and non-equilibrium 
entropy function was a regular function of the temperature; nevertheless that 
model involved a temperature dependence of the volume jump near critical point 
in the form Ai' ^ \6 — 6c\. Furthermore, the latent heat does not approach to 
zero at critical point. 

This shortcoming can be avoided in ad hoc way assuming a cancelling diver- 
gence contribution in rjo or following more closely the Landau theory of second 
order phase transitions, in which the entropy remains a regular function at the 
critical point, while the equilibrium value of the order parameter has a square 
root singularity at critical point. 
We suggest to consider the potential 

f{p, 9, ip) = /o(p, 9) - a{ui9) + 1) ln(l - ip^) - a^"" - 2h{p, 9)^, (59) 

with a > and u{9) an increasing function of 9 such that 

u(6l) e (-l,+oo), u(6lc) = 0, u(0+) = -l, (60) 

as, for example, in the following class of functions: 

91 2^ 



u{9) 



1- ^ 



g > 0, /? > 0. (61) 



As before, h has to satisfy h{po{9), 9) — and ft. ^ for p ^ Pq{S). 
The minima are solutions of the third order algebraic equation 

if^ +uip=-{l~-ip^). (62) 
a 



11 



The same algebraic equation for the minima would be obtained (except for 
solutions (p > ±1) from the fourth order polynomial potential 



if^ h f Lp'^ 



(63) 



which is the classical form of the Ginzburg-Landau potential for second order 
phase transitions. The employment of the logarithmic potential (also used in 
different contexts, see for example [9 ) has the advantage that the only attainable 
minima are constrained in the interval (—1,1). This would not be of great 
importance for the only study of critical region, where order parameter is small: 
</? <C 1. On the contrary, if one wants the phase field to be bounded under all 
the temperature and pressure conditions, the fourth order polynomial potential 
would not be suitable. If m < 0, there are two local minima — 1 < i^i < < 

h/a 





-a 



Figure 3: Structure of the minima of the potential f{p,8,ip) at varying h > 0,u > 0. In the 
region h < the structure is the same with the interchange of ip with —ip. 

(^2 < 1 for \h/a\ < h{u), while there is only one local minimum —1 < (p < 1 ii 
\h/a\ > h{u). (The function h{u) has not a simple expression, but 2|u/3|'^''^ < 
h < \u\; see Fig. [3]). For u > 0, there is always only one minimum in [—1, 1]. 
In particular, consider the case h =^ 0; ior 6 < 9c (on the coexistence line), the 
minima are given by 

^±=±Hi/2e[-l,l], (64) 

while for 9 > 9c (on the critical isochore) the only minimum is = 0. 

For ft. 7^ 0, the minima have complicated algebraic expressions which are not very 
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illuminating; we only remark that the absolute minimum satisfies the following 
estimate: 

I^ Lfr'^l , <l~^^< i,^"'""' for~l<u<0. (65) 
' ' |ft|/a+ |u| - ^ - \h\/a + l' ^ ' 

In Fig. [4] it is plotted the value of (p in function of h for three values of u. The 




Figure 4: Equilibrium solution as a function of h/a for three values of n: u = —0.4 (dashed 
line), li = (continuous line), u = 0.4 (dash-dotted line). For u < and \h/a\ < h there are 
two minima, which have equal height at /i = 0; this gives rise to the characteristic hysteresis 
loop of the dashed line. 

equation for volume is 

V = (/o)p - 2hpLp, (66) 
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which, owing to algebraic equation (62) for the equilibrium value implies the 
following state equation 

(/op - + ^(/op - + ^uhliU, - ^) = (67) 

The specific volume jump at the transition is 

Av^Ahp\u\^'^. (68) 

If u ~ {9/6c~l), the quantity Av has the typical mean-field singularity |6'— 0^^^^ 
near the critical point, needlessly to have u{0) also singular. 
Finally, equilibrium entropy is given by 

V = -ih)e + au'{e) ln(l - ^'') + 2he^. (69) 

Now, u{6) is a regular function of f?, and so f] is not diverging in the neighbour- 
hood of the critical point. 

With this choice of the Ginzburg-Landau potential, the equilibrium value (ab- 
solute minimum) of the order parameter (p is not constant in the inner of the 
regions of the phase diagram marked as 'liquid' and 'vapour'. In this sense, 
the Ginzburg-Landau equation is never completely decoupled from the Navier- 
Stokes equation. However, as 

l-(?2<^(l-|u|), (70) 

the order parameter is approximatively ±1 in the region \h\ 3> a (far from 
transition line) and/or |it| ^ 1 (well below the critical temperature). 
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